home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-04
/
amsf20.zip
/
AMST3.FOR
< prev
next >
Wrap
Text File
|
1992-01-06
|
7KB
|
233 lines
C ******************************************************************
C * *
C * E I G E N *
C * *
C ******************************************************************
PROGRAM EIGEN
IMPLICIT INTEGER*4(I-N)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON MAVAIL,IA(30000)
MAVAIL=30000
C ... SAMPLE: EIGEN PROBLEM SOLVER
PRINT *,' '
PRINT *,'EIGEN PROBLEM SOLVER'
PRINT *,' '
PRINT *,'ENTER EQUATION NUMBER '
READ *, N
PRINT *,'NOW ENTER STIFFNESS MATRIX & MASS MATRIX'
PRINT *,' '
C ... TEST IN-CORE DATA MANAGEMENT
CALL DBOPEN(1,'TSTDAT.DAT','NEW')
PRINT *,'DATABASE OPENED.'
CALL DEFINE(1,'STIF',0,1,N,N,0,NA)
CALL DEFINE(1,'MASS',0,1,N,N,0,NB)
CALL DEFINE(1,'VCTR',0,1,N,N,0,NC)
CALL DEFINE(1,'EIGN',0,1,N,1,0,ND)
CALL DEFINE(1,'WORK',0,1,N,1,0,NE)
print *,'DEFINE MATRICES DONE.'
CALL MATINP(1,'STIF')
CALL MATINP(1,'MASS')
CALL MATOUT(1,'STIF')
CALL MATOUT(1,'MASS')
RTOL = 1.0D-12
CALL JACOBI(IA(NA),IA(NB),IA(NC),IA(ND),IA(NE),N,RTOL,15,0,6)
CALL MATOUT(1,'EIGN')
CALL MATOUT(1,'VCTR')
CALL DIR(0)
CALL DBCLOS(1,'DELETE')
STOP 'DONE.'
END
C ******************************************************************
C * *
C * J A C O B I *
C * *
C ******************************************************************
SUBROUTINE JACOBI (A,B,X,EIGV,D,N,RTOL,NSMAX,IFPR,IOUT)
IMPLICIT INTEGER*4(I-N)
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A(N,N),B(N,N),X(N,N),EIGV(N),D(N)
C
C TO SOLVE THE GENERALIZED EIGENPROBLEM USING THE GENERALIZED
C JACOBI ITERATION
C
C INPUT VARIABLES :
C A(N,N) = STIFFNESS MATRIX (ASSUME POSITIVE DEFINITE)
C B(N,N) = MASS MATRIX (ASSUME POSITIVE DEFINITE)
C D(N) = WORKING VAECTOR
C N = ORDER OF MATRIX A AND B
C RTOL = CONVERGENCE TOLERANCE (SET TO 1.E-12)
C NSMAX = MAXIMUM NUMBER OF SWEEP ALLOWED (SET TO 15)
C IFPR = FLAG FOR PRINTING DURING ITERATION
C ( EQ.0 NO PRINTING, EQ.1 PRINT INTERMEDIATE RESULTS)
C IOUT = OUTPUT LUN
C OUTPUT VARIABLES :
C A(N,N) = DIAGONALIZED STIFFNESS MATRIX
C B(N,N) = DIAGONALIZED MASS MATRIX
C X(N,N) = EIGENVECTORS STORED COLUMNWISE
C EIGV(N) = EIGENVALUES
C
C INITIALIZE EIGENVALUE AND EIGENVECTOR MATRICES
DO 10 I=1,N
IF (A(I,I).LE.0..OR.B(I,I).LT.0.) THEN
WRITE(IOUT,2020)
STOP 'JACOBI'
END IF
D(I)=A(I,I)/B(I,I)
10 EIGV(I)=D(I)
DO 30 I=1,N
DO 20 J=1,N
20 X(I,J)=0.
30 X(I,I)=1.0
IF(N.EQ.1) RETURN
C INITIALIZE SWEEP COUNTER AND EIGEN ITERATION
NSWEEP=0
NR=N-1
C
C WE START ITERATION
40 NSWEEP=NSWEEP+1
IF(IFPR.EQ.1) WRITE(IOUT,1000) NSWEEP
C CHECK IF PRESENT OFF-DIAGONAL ELEMENT IS LARGE ENOUGH TO
C REQUIRE ZEROING
EPS=(0.01**NSWEEP)**2
DO 210 J=1,NR
JJ=J+1
DO 210 K=JJ,N
TT=A(J,K)*A(J,K)
TB=A(J,J)*A(K,K)
EPTOLA=DABS(TT/TB)
TT=B(J,K)*B(J,K)
TB=B(J,J)*B(K,K)
EPTOLB=TT/TB
IF ((EPTOLA.LT.EPS).AND.(EPTOLB.LT.EPS)) GO TO 210
C IF ZEROING IS REQUIRED, CALCULATE THE ROTATION MATRIX ELEMENT CA, CG
AKK=A(K,K)*B(J,K)-B(K,K)*A(J,K)
AJJ=A(J,J)*B(J,K)-B(J,J)*A(J,K)
AB=A(J,J)*B(K,K)-A(K,K)*B(J,J)
CHECK=(AB*AB+4.0*AKK*AJJ)/4.0
IF (CHECK) 60,60,70
60 WRITE (IOUT,1004) CHECK
STOP
70 SQCH=DSQRT(CHECK)
D1=AB/2.0+SQCH
D2=AB/2.0-SQCH
DEN=D1
IF (DABS(D2).GT.DABS(D1)) DEN=D2
IF (DEN) 90,80,90
80 CA=0.
CG=-A(J,K)/A(K,K)
GO TO 100
90 CA=AKK/DEN
CG=-AJJ/DEN
C
C WE PERFORM THE GENERALIZED ROTATION
100 IF (N-2) 101,180,101
101 JP1=J+1
JM1=J-1
KP1=K+1
KM1=K-1
C
IF (JM1-1) 120,110,110
110 DO 105 I=1,JM1
AJ=A(I,J)
BJ=B(I,J)
AK=A(I,K)
BK=B(I,K)
A(I,J)=AJ+CG*AK
B(I,J)=BJ+CG*BK
A(I,K)=AK+CA*AJ
105 B(I,K)=BK+CA*BJ
C
120 IF (KP1-N) 130,130,140
130 DO 125 I=KP1,N
AJ=A(J,I)
BJ=B(J,I)
AK=A(K,I)
BK=B(K,I)
A(J,I)=AJ+CG*AK
B(J,I)=BJ+CG*BK
A(K,I)=AK+CA*AJ
125 B(K,I)=BK+CA*BJ
C
140 IF (JP1-KM1) 150,150,180
150 DO 160 I=JP1,KM1
AJ=A(J,I)
BJ=B(J,I)
AK=A(I,K)
BK=B(I,K)
A(J,I)=AJ+CG*AK
B(J,I)=BJ+CG*BK
A(I,K)=AK+CA*AJ
160 B(I,K)=BK+CA*BJ
180 AK=A(K,K)
BK=B(K,K)
A(K,K)=AK+2*CA*A(J,K)+CA*CA*A(J,J)
B(K,K)=BK+2*CA*B(J,K)+CA*CA*B(J,J)
A(J,J)=A(J,J)+2*CG*A(J,K)+CG*CG*AK
B(J,J)=B(J,J)+2*CG*B(J,K)+CG*CG*BK
A(J,K)=0.0
B(J,K)=0.0
C
C UPDATE EIGENVECTORS
DO 190 I=1,N
XJ=X(I,J)
XK=X(I,K)
X(I,J)=XJ+CG*XK
190 X(I,K)=XK+CA*XJ
C
210 CONTINUE
C
DO 220 I=1,N
220 EIGV(I)=A(I,I)/B(I,I)
IF(IFPR.EQ.0) GO TO 227
WRITE (IOUT,1005)
WRITE (IOUT,1002) (EIGV(I),I=1,N)
227 CONTINUE
C
C CHECK FOR CONVERGENCE
DO 240 I=1,N
TOL=RTOL*D(I)
DIF=DABS(EIGV(I)-D(I))
IF(DIF.GT.TOL) GO TO 300
240 CONTINUE
C
C CHECK IF ALL OFF-DIAG ELEMENTS ARE SATISFACTORILY SMALL
EPS=RTOL**2
DO 260 J=1,NR
JJ=J+1
DO 260 K=JJ,N
TT=A(J,K)*A(J,K)
TB=A(J,J)*A(K,K)
EPSA=DABS(TT/TB)
TT=B(J,K)*B(J,K)
TB=B(J,J)*B(K,K)
EPSB=TT/TB
IF ((EPSA.LT.EPS).AND.(EPSB.LT.EPS)) GO TO 260
GO TO 300
260 CONTINUE
C
DO 310 I=1,N
DO 310 J=I,N
B(J,I)=B(I,J)
310 A(J,I)=A(I,J)
RETURN
C
300 DO 320 I=1,N
320 D(I)=EIGV(I)
IF (NSWEEP.LT.NSMAX) GO TO 40
DO 330 I=1,N
DO 330 J=I,N
B(J,I)=B(I,J)
330 A(J,I)=A(I,J)
RETURN
C
1000 FORMAT (27H0SWEEP NUMBER IN *JACOBI* =, I4)
1002 FORMAT (1H ,12E11.4)
1004 FORMAT (37H0***ERROR SOLUTION STOP IN *JACOBI*, / 12X,
1 8HCHECK = , E20.14 / 1X)
1005 FORMAT (36H0CURRENT EIGENVALUES IN *JACOBI* ARE, / 1X)
2020 FORMAT ('0**ERROR** MATRIX NOT POSITIVE DEFINITE IN JACOBI'/)
C
END